#include <iostream>
#include <string>
#include <vector>
#include <fstream>

// include libigl headers

#include <igl/grad.h>
#include <igl/jet.h>
#include <igl/per_vertex_normals.h>
#include <igl/opengl/glfw/Viewer.h>
#include <igl/principal_curvature.h>
#include <igl/jet.h>

// include slTools headers
#include "particle2realflowUtil/slVector.H"
#include "particle2realflowUtil/kdTree.H"


// include my private functions
#include "myFunctions.h"

using namespace Eigen;
using namespace std;

Eigen::MatrixXd V, U, V_t, V_d, C;
Eigen::MatrixXi F;
Eigen::SparseMatrix<double> L;
Eigen::SparseMatrix<double> Weight_L, Weight_V;
Eigen::VectorXd GussianCurvature;
std::vector<SlVector3> Particles;
Eigen::MatrixXd ParticleAnisotropy;
Eigen::MatrixXd VerticesAnisotropy;
Eigen::MatrixXd VerticesOffset;
Eigen::VectorXi numPointsOfVer;
KDTree * kdTreePointer = NULL;

arguConfigure argumentInfo;
bool isVisualize = true;

int main(int argc, char *argv[])
{
	igl::opengl::glfw::Viewer viewer;
	float inwards = 0.0;

	readConfigureFile("D:/Workspace/Graphics/libigl-example-project/configure.txt", argumentInfo);

	if (argumentInfo.fileType == 1) {
		SkinningParticles( (argumentInfo.workDirectory+argumentInfo.fileName+".txt").c_str(), inwards, V, F, Particles, argumentInfo);
	}
	else if (argumentInfo.fileType == 2 ) {
		readModelFile(argumentInfo.workDirectory + argumentInfo.fileName + ".obj", V, F);
		readParticlData( argumentInfo.workDirectory + argumentInfo.fileName + ".txt", Particles);
	}
	else if (argumentInfo.fileType == 3) {
		readModelFile(argumentInfo.workDirectory + argumentInfo.fileName + ".off", V, F);
		readParticlData(argumentInfo.workDirectory + argumentInfo.fileName + ".txt", Particles);
	}
	else {
		std::cout << "read configure file failed." << endl;
		return 0;
	}


	outputModelinfo(V, F);

	//estimate the size of the model
	outputModelLength(V);

	// Compute centroid and subtract, scale it to uniformed scale
	// but now no need to do this.

	//estimate the size of the uniformed model

	assert(V.cols()==3);

	V_t = V;



	/*
	// Alternative construction of same Laplacian
	SparseMatrix<double> G, K;
	// Gradient/Divergence
	igl::grad(V, F, G);
	// Diagonal per-triangle "mass matrix"
	VectorXd dblA;
	igl::doublearea(V, F, dblA);
	// Place areas along diagonal #dim times
	const auto & T = 1.*(dblA.replicate(3, 1)*0.5).asDiagonal();
	// Laplacian K built as discrete divergence of gradient or equivalently
	// discrete Dirichelet energy Hessian
	K = -G.transpose() * T * G;
	cout << "|K-L|: " << (K - L).norm() << endl;
	*/


	const auto &key_down = [](igl::opengl::glfw::Viewer &viewer, unsigned char key, int mod)->bool
	{
		switch (key)
		{
		case 'r':
		case 'R':
			viewer.data().clear();
			U = V;
			viewer.data().set_mesh(U, F);
			break;
		case 'w':
		case 'W': 
		{
			writeIntoOBJ("D:/Workspace/Graphics/model/newModel.obj", U, F);
			break; 
		}

		case 'i':
		case 'I':
		{
			viewer.data().clear();
			viewer.data().set_mesh(U, F);
			break;
		}

		case 'k':	// build KD-tree and compute vertices singularities
		case 'K':
		{
			if (kdTreePointer != NULL) {
				delete kdTreePointer;
			}

			cout << "Number of particles: " << Particles.size() << endl;

			//establish kd-tree
			kdTreePointer = new KDTree(Particles);
			cout << "KD-Tree has been established!" << endl;

			ParticleAnisotropy.resize(Particles.size(), 3);

			calculateParticleAnisotropy(Particles, *kdTreePointer, argumentInfo.pcaMaxNumParticleNeighbors , 
				argumentInfo.h*0.866* argumentInfo.pcaParticleRadiusRate, ParticleAnisotropy);


			VerticesAnisotropy.resize(U.rows(),1);
			numPointsOfVer.resize(U.rows());

			calculateMeshAnisotropy(U, Particles, *kdTreePointer, argumentInfo.pcaMaxNumVerticesNeighbors,
				 argumentInfo.h*0.866 * argumentInfo.pcaVerticesRadiusRate , ParticleAnisotropy, 
				VerticesAnisotropy, numPointsOfVer);


			break;
		}
		case 'P':
		case 'p':	// principal component analysis
		{
			if (kdTreePointer == NULL) {
				std::cout << "Error! KD-Tree has not been established!" << endl;
				break;
			}

			ParticleAnisotropy.resize(Particles.size(), 3);

			calculateParticleAnisotropy(Particles, *kdTreePointer, argumentInfo.pcaMaxNumParticleNeighbors,
											argumentInfo.h*0.866 * argumentInfo.pcaParticleRadiusRate, ParticleAnisotropy);


			VerticesAnisotropy.resize(U.rows(), 1);

			calculateMeshAnisotropy(U, Particles, *kdTreePointer, argumentInfo.pcaMaxNumVerticesNeighbors,
					argumentInfo.h*0.866 * argumentInfo.pcaVerticesRadiusRate , 
					ParticleAnisotropy, VerticesAnisotropy, numPointsOfVer);


			break;
		}
			
		case ' ':	//minimize quadratic energy
		{

			{
				// Compute Laplace-Beltrami operator: #V by #V
				getLaplacianMatrix(U, F, L);

				// Solve (M-delta*L) U' = M*U
				//const auto & S = (I - 0.0002*L);
				SparseMatrix<double, RowMajor> S(2 * L.rows(), L.cols());
				{
					// weight multiply the Identity Matrix
					SparseMatrix<double, RowMajor> I(L.rows(), L.cols());
					I.setIdentity();

					SparseMatrix<double, RowMajor> LaplacianWeightMatrix(L.rows(), L.cols());
					cout << "Laplacian Matrix size: " << L.size() << endl;
					cout << "Laplacian weight Matrix rows: " << LaplacianWeightMatrix.rows() << endl;
					cout << "Laplacian weight Matrix cols: " << LaplacianWeightMatrix.cols() << endl;

					getWeightMatrix(VerticesAnisotropy.col(0), argumentInfo.LaplacianWeight,
						argumentInfo.LaplacianWeightRate , LaplacianWeightMatrix);


					S.topRows(L.rows()) = (LaplacianWeightMatrix*L).topRows(L.rows());
					S.bottomRows(L.rows()) = argumentInfo.VerticesWeight*I.bottomRows(L.rows());
				}
				cout << "Matrix S: " << S.rows() << " * " << S.cols() << endl;
				MatrixXd Beta(2 * L.rows(), 3);
				Beta << MatrixXd::Zero(L.rows(), 3), argumentInfo.VerticesWeight * U;
				cout << "Matrix Beta: " << Beta.rows() << " * " << Beta.cols() << endl;

				solveLinearLeastEquation(S, Beta, U);
			}
			
			// output the size of the mesh after deformed
			outputModelLength(U);
			

			break;
		}
		case '1': // Compute vertices offset through particle layout
		{
			if (kdTreePointer == NULL) {
				std::cout << "Error! KD-Tree has not been established!" << endl;
				break;
			}
			

			break;
		}
		case '2': // Compute vertices MCF wieght Weight_L & Weight_V
		{
			if (kdTreePointer == NULL) {
				std::cout << "Error! KD-Tree has not been established!" << endl;
				break;
			}

			break;
		}
		case '4':
		{
			V_d = U;

			viewer.data().add_edges(V_d, V_t, Eigen::RowVector3d(1.0, 0.0, 0.0));

			Eigen::MatrixXd delta_result = V_d - V_t;
			Eigen::VectorXd errorVec = delta_result.rowwise().norm();
			double meanError = errorVec.mean();
			double maxError = errorVec.maxCoeff();
			cout << "The mean error is: " << meanError << endl;
			cout << "The max error is: " << maxError << endl;

			break;
		}
		case '6': // use mean curvature to dye vertices
		{
			Eigen::VectorXd LaplacianMeanCurvature = 0.5*(-L*U).rowwise().norm();

			double averageMC = LaplacianMeanCurvature.sum()/ LaplacianMeanCurvature.size();
			cout << "Laplace:Average mean curvature is: " << averageMC << endl;
			cout << "Laplace:Maximum mean curvature is: " << LaplacianMeanCurvature.maxCoeff() << endl;
			cout << "Laplace:Minimum mean curvature is: " << LaplacianMeanCurvature.minCoeff() << endl;

			
			igl::jet( (LaplacianMeanCurvature + VectorXd(LaplacianMeanCurvature.size()).setConstant(1.0))
									.array().log().matrix()  , true, C);
			viewer.data().set_colors(C);
			break;
		}
		case '7': // use mean curvature to dye vertices
		{
			MatrixXd PD1, PD2;
			VectorXd PV1, PV2;
			igl::principal_curvature(U, F, PD1, PD2, PV1, PV2);
			Eigen::VectorXd meanCurvature = (PV1+PV2)*0.5;

			double averageMC = meanCurvature.sum() / meanCurvature.size();
			cout << "Laplace:Average mean curvature is: " << averageMC << endl;
			cout << "Maximum mean curvature is: "<< meanCurvature.maxCoeff() <<endl;
			cout << "Minimum mean curvature is: " << meanCurvature.minCoeff() << endl;
			igl::jet(meanCurvature, true, C);
			viewer.data().set_colors(C);
			break;
		}
		case '8': // use the maximum difference of SVD singularity value to dye vertices
		{
			if (kdTreePointer == NULL) {
				std::cout << "Error! KD-Tree has not been established!" << endl;
				break;
			}

			Eigen::VectorXd deltaSingularity = VerticesAnisotropy.col(0);
			igl::jet(deltaSingularity, true, C);
			viewer.data().set_colors(C);
			break;
		}
		case '9': // use the num of points of per vertices to dye
		{
			if (kdTreePointer == NULL) {
				std::cout << "Error! KD-Tree has not been established!" << endl;
				break;
			}

			Eigen::VectorXd deltaSingularity = VerticesAnisotropy.col(0);
			igl::jet(numPointsOfVer, true, C);
			viewer.data().set_colors(C);
			break;
		}
		case 'c': 
		case 'C':	// control: 
		{
			//clean the buffer
			

			int controlNum = 0;
			cout << "please input control number:" << endl;
			cout << "1 represents for searching area:" << endl;
			cin >> controlNum;
			if (controlNum == 1) {
				cout << "please input PCA particle neighbor area:" << endl;
				cin >> argumentInfo.pcaParticleRadiusRate;
				cout << "please input PCA vertices neighbor area:" << endl;
				cin >> argumentInfo.pcaVerticesRadiusRate;
			}
			else {
			}

			cout << "control finished." << endl;
			break;
		}
		default:
			return false;
		}
		// Send new positions, update normals, recenter
		viewer.data().set_vertices(U);
		//viewer.data().compute_normals();
		viewer.core.align_camera_center(U, F);
		return true;
	};


	// Initialize smoothing with base mesh
	U = V;
	viewer.data().set_mesh(U, F);
	viewer.data().show_lines = false;
	getLaplacianMatrix(U, F, L);
	//viewer.data().set_colors(C);
	viewer.callback_key_down = key_down;

	cout << "Press [space] to smooth." << endl;
	cout << "Press [r] to reset." << endl;
	return viewer.launch();
}